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ABSTRACT 

An explanation for the large random velocities of pulsars is presented. Like many other models, we propose 
that the momentum imparted to the star is given at birth. The ultimate source of energy is provided by the intense 
optically thick neutrino flux that is responsible for radiating the proto-neutron star's gravitational binding energy 
during the Kelvin-Helmholtz phase. The central feature of the kick mechanism is a radiative-driven magnetoa- 
coustic instability, which we refer to as "neutrino bubbles." Identical in nature to the photon bubble instability, the 
neutrino bubble instability requires the presence of an equilibrium radiative flux as well as a coherent steady back- 
ground magnetic field. Over regions of large magnetic flux densities, the neutrino bubble instability is allowed to 
grow on dynamical timescales ^ 1ms, potentially leading to large luminosity enhancements and density fluctua- 
^ ' tions. Local luminosity enhancements, which preferentially occur over regions of strong magnetic field, lead to a 

Q I net global asymmetry in the neutrino emission and the young neutron star is propelled in the direction opposite to 

these regions. For favorable values of magnetic field structure, size, and strength as well as neutrino bubble satu- 
\^ ' ration amplitude, momentum kicks in excess of 1000 kms"' can be achieved. Since the neutrino-powered kick is 

delivered over the duration of the Kelvin-Helmholtz time ^ a few seconds, one expects spin-kick alignment from 
this neutrino bubble powered model. 

^ , Subject headings: stars:evolution-stars:oscillations-MHD-instabilities- gravitational waves-pulsars:general 

1. INTRODUCTION 

. Shortly after their discovery, it was understood that neutron stars travel with random velocities significantly in excess of the 
massive stars from which they were born (Gunn & Ostriker 1971). More recently, studies of pulsar velocity distributions suggest that 
\ neutron stars indeed possess large proper motions and that the momentum kick is most likely given at birth possibly with a bimodal 
velocity distribution (Lyne & Lorimer 1994; Lorimer, Bailes, & Harrison 1997; Fryer, Burrows, & Benz 1998; Cordes & Chernoff 
1998; Arzoumanian, Chernoff, & Cordes 2002). The orbital behavior of individual double neutron star and black hole systems also 
imply that neutron stars receive natal kicks (Fryer & Kalogera 1997; Kramer 1998; Tauris et al. 1999; Wex, Kalogera, & Kramer 
Q 2000; Mirabel et al. 2002). 

i—{ From a theoretical point of view, a natal kick is a reasonable outcome of stellar collapse and the subsequent explosion. Approxi- 
c/2 , mately 10^' erg of kinetic energy is released by the prompt explosion mechanism, which can be thought of as as shell of ~ 10~'Mo 
Ci ' traveling radially outward at ~ 10"* cms"'. Thus, a momentum asymmetry in the shock of ^ 10% could potentially endow a neutron 
star with a kick velocity ~ 10^ cms"' = 10^ kms"'. Motivated by this simple kinematic argument. Burrows & Hayes (1996) pursued 
a "hydrodynamically-driven" kick mechanism by artificially deforming the outer portions of an Fe core progenitor and numerically 
/\ ' evolving the subsequent collapse and bounce. They found that a "mass-rocket" produced a kick of ~ 500 kms"'. One possibility for 
H a substantial seed asymmetry may be due to unstable g-modes in progenitor Fe core, which are driven by the e-mechanism in the 
. F^. . outer shell (Lai & Goldreich 2000). 

The intense optically thick neutrino radiation is another potential reservoir of energy for natal kicks. During the Kelvin-Helmholtz 
phase, neutrino emission is responsible for radiating ~ 99% of the star's gravitational binding energy, which amounts to ~ a few 
X 10^^ erg. Asymmetric neutrino surface emission of order ~ a few% leads to a kick of ~ lO^'kms"'. The magnetic field of a 
proto-neutron star may provide the needed systematic asymmetry in such a profoundly spherically symmetric environment, in which 
the neutrinos can couple. Since the weak interaction violates parity conservation, the neutrinos preferentially scatter towards one 
magnetic pole over the other. This "neutrino/magnetic-driven" mechanism requires global magnetic field strengths ^ 10"'G in order 
to produce a kick of - 10^ kms"' (Arras & Lai 1999a, 1999b). 

In this paper, we describe an alternative model for neutron star kicks that can be viewed, in the broadest sense, as a synthesis of the 
hydrodynamic- and neutrino/magnetic- driven models. By considering general radiative envelopes that are stratified, optically thick, 
magnetized, and highly electrically conducting, Blaes & Socrates (2003, hereafter BS03) showed that local magnetoacoustic pertur- 
bations are susceptible to radiative driving, a phenomenon known as "photon bubbles" (Arons 1992; Gammie 1998). The non-linear 
evolution of the photon bubble instability is far from certain as research of this phenomenon is in the developmental stages (Begelman 
2001, Turner et al. 2004). Regardless, it is reasonable to suspect that radiative flux enhancements occur in the saturated state since 
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the radiative flux provides the free energy for the instability. In the case of proto-neutron star atmospheres, neutrino bubbles rather 
than photon bubbles may vigorously operate in surface regions that are substantially magnetized since neutrinos rather than photons 
are the particle species responsible for the diffusive transfer of energy. Over these regions, neutrino flux enhancements propel the 
star in the opposite direction, potentially providing a natal kick (See Figure 1). In a way, one can view this mechanism as a "neutrino 
magnetohydrodynamically-driven" mechanism. The main focus of this work is to understand the mechanics, thermodynamics, and 
stability properties of local radiation magnetohydrodynamic oscillations in proto-neutron stars so that we can assess the feasibility, 
to first approximation, of a neutrino bubble powered kick mechanism. 

The organization of this paper is as follows. In §2 we provide some initial theoretical arguments in order to motivate this work. The 
importance of radiation-driven fluid instabilities is stressed by comparing the envelope structure of proto-neutron stars with those on 
the main sequence. Based on the magnetic field structure of the solar convection zone, we describe magnetic field configurations on 
the surface of proto-neutron stars that are favorable for neutrino bubble powered kicks. A detailed analysis of local magnetoacoustic 
perturbations in the limit of rapid neutrino diffusion is carried out in §§3-5. Our simplifying assumptions and problem set up are 
sketched out in §3. The thermodynamic behavior of short-wavelength oscillations in degenerate radiating fluids is investigated in §4. 
The driving mechanism of the neutrino bubble instability is presented §5 in terms of the relationship between pressure and density 
as well as force and velocity. Also in §5, the neutrino bubble stability criterion is evaluated in terms of parameters relevant for proto- 
neutron star surfaces. A Monte Carlo radiation transfer study is performed in §6 where we estimate the value of the kick as well as 
the gravitational wave strain produced by the globally asymmetric neutrino radiation. A comparison to other kick mechanisms and a 
critical summary are given in §7. 

2. THE KICK MECHANISM 

A proto-neutron star (PNS) releases ^ 99% of the accumulated binding energy from its Fe core progenitor in the form of neutrinos. 
The neutrinos escape the PNS nearly uniformly across all neutrino flavors during the Kelvin-Helmholtz phase, which lasts for ~ 10 s 
after its birth. If the neutrino emission is asymmetric by roughly a few percent, the PNS receives a natal kick ^ 10"^kms~' (Arras & 
Lai 1999a). Since the neutrino emission by far dominates the energetics during the Kelvin-Helmholtz phase, the main difficulty in 
formulating a neutrino-powered kick model rests in identifying an intrinsic systematic asymmetry in which the neutrinos can couple. 

A natural candidate for inducing a neutrino asymmetry is the magnetic field of the PNS. On their own, magnetic fields alter 
neutrino-nucleon scattering and absorption via coupling between the background field and the nucleon's magnetic dipole moment. 
Since the weak interaction violates parity, the introduction of an axial magnetic field into the scattering matrix causes the neutrinos 
to preferentially scatter toward one magnetic pole over the other. However, this effect exists only as long as neutrinos are allowed to 
deviate from local thermal and chemical equilibrium (Arras & Lai 1999a). At and above this layer, neutrinos begin to mechanically 
decouple from the matter as well, implying a lack of momentum exchange between the neutrinos and the entire PNS. Thus, strong 
~ 10'^ G global magnetic fields are required to induce a significant neutrino asymmetry above the energy decoupling layer. 

Here, we also exploit the magnetic field for the purpose of generating neutrino-powered natal kicks. Rather than altering the 
thermal transport properties of neutrinos, we propose that the magnetic field induces a radiation-driven fluid instability at and below 
the neutrinosphere. The fluid instability is acoustic in nature and is driven in the exact same way as the "photon bubble" instabilities, 
which may operate in neutron star accretion caps, relativistic accretion disks, and magnetized main sequence stellar envelopes (Arons 
1992, Gammie 1998, BS03). In the case of PNSs, neutrinos are the species with relatively large mean free paths and therefore enable 
the diffusive transport of energy. It follows that in the envelopes of PNSs, the neutrino bubble instability is the relevant mechanism 
when considering radiative magnetosonic stability. As the instability grows and saturates, a luminosity perturbation is induced about 
regions of large magnetic flux density and the PNS is propelled in the direction opposite to these "starspots." 

In order for a radiation-driven fluid instability to saturate with a large luminosity perturbation, one expects the radiation pressure 
to be at least comparable in value to the gas pressure. For PNSs, the average neutrino luminosity during the Kelvin-Helmholtz era is 
~ 0.01 the Eddington value, implying that the equilibrium radiation pressure is ~ 0.01 the equilibrium gas pressure. However, at and 
below the neutrinosphere, the radiation and neutrinos are in thermal and chemical equilibrium with one another on the dynamical 
timescales of interest, which are of order 1 ms - the acoustic crossing time over a gas pressure scale height. We will demonstrate 
that the spatial distribution of both temperature and neutrino chemical potential for a given acoustic MHD mode are responsible 
for driving the neutrino bubble instability. Since the the temperature and relative composition of the gas are locked to that of the 
neutrinos, the gas itself behaves in a fashion similar to the radiation and provides for the majority of the driving. It follows that at 
least to linear order, the stress responsible for radiative driving is comparable in magnitude to the acoustic restoring force, opening 
up the possibility for large saturation amplitudes. 

Radiation-driven acoustic instabilities, e.g. the K-mechanism, are ultimately driven by the background radiation flux coupling to 
the fluid displacement and density perturbation. Yet, damping from radiative diffusion only depends on the total pressure of the gas. 
The larger the ratio of flux to pressure, the more likely driving will overcome damping. In other words, as the radiation field becomes 
increasingly anisotropic, the free energy per quanta available for acoustic driving increases accordingly. Due to their stiff equation 
of state above nuclear saturation density and the relatively low number of total optical depths ~ a few thousand, the radiation fields 
of PNSs are highly anisotropic in comparison to those of main sequence stars whose luminosities are ~ 0.01 Ledd in photons. To 
illustrate this point, compare the fraction of stellar mass lying one optical depth below the surface of last scattering for a PNS and a 
4Mq A-star 

™-™^ -> ^ s ^ -"-^ m m iw^) 



3 



V 




Fig. 1 . — Sketch of the kick mechanism. Starspot-like magnetic field structures are generated at the base of convection zone and are transported to the upper-base 
of the convection zone by both buoyancy and turbulent pumping in a few overturn times. If the presence of strong, long-lived, large scale magnetic "starspots" can 
induce neutrino bubble instability for a significant fraction of the Kelvin-Helmholtz time, the proto-neutron star will receive a natal kick due to enhanced neutrino 
transport in the vicinity of the starspot. For average initial rotation rates and random starspot placement, the rotation axis will most likely be aligned with the kick 
for a given neutron star. 



where p,-, and is the radius of the surface of last scattering, the scale height of the surface of last scattering, the density at 

the last scattering surface, and the central core density, respectively. From the above argument, one suspects that a greater fraction 
of the total stellar mass is responsible for radiative-acoustic driving in PNSs than in main sequence stars of comparable Eddington 
factors. Thus, PNSs are promising sights for radiation-driven acoustic instability of any type since a larger mass fraction responsible 
for driving implies a relatively large luminosity perturbation. If, for example, the K-mechanism (Baker & Kippenhanh 1962, Unno et 
al. 1989) or strange modes (Glatzel 1994) operate in the envelopes of PNSs, large luminosity perturbations may follow. An increase 
in the overall integrated luminosity of ~ 10% during the first few hundred milliseconds may significantly contribute to the explosion 
mechanism. Since the focus of this work is centered on producing a potential kick mechanism, our analysis will be heavily weighted 
toward radiation-MHD driving mechanisms which are only located on patches of large magnetic flux density. 

The nature of the magnetic field structures in the vicinity of the neutrinosphere in part, determines whether or not neutrino bubble 
instability can power kicks, or even explosions. In the following discussion, mainly phenomenological arguments are put forth in 
order to motivate the possibility that coherent magnetic formations can exist for a substantial fraction of the Kelvin-Helmholtz time 
while covering a sizable fraction of the star, with large (10'"*- lO'^G) field strengths. 



2.1. Convection and Magnetism in Proto-Neutron Stars 
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It is believed that PNSs convect as they cool. At very early times the escaping shock dissociates ambient nuclei into baryons, 
heating the mantle while producing a negative entropy gradient that drives convection (Bethe, Brown & Cooperstein 1987; Wilson 
& Mayle 1988). Due to their large energies and mean free paths, neutrinos radiate the entropy excess left over by the prompt shock, 
immediately quenching the convection in < 100 ms. Soon after this somewhat superficial episode of convective overturn, a PNS 
suffers from deep convective instability for the duration of the Kelvin-Helmholtz time (Thompson & Murray 2001, hereafter TMOl; 
Burrows & Lattimer 1986 Burrows 1987; TD93; Keil, Janka, & Mueller 1996, hereafter KJM96; Pons et al. 1999). KJM96 find that 
Ledoux-type convection quickly recedes from the neutrinosphere in ~ 20-30 ms, moving deep into the core and at about ^ Is after 
bounce, the entire star is convective. TD93 point out that in PNSs, the ratio of convective kinetic energy to gravitational binding 
energy ~ 100 times larger than its progenitor during the Silicon burning phase and ^ 10^- lO'* times larger when the progenitor 
resided on the main sequence. In fact, when compared this way, convection in PNSs is ^ 10"* times more violent than the Sun. 
Primarily motivated by this fact, TD93 conclude that PNSs are efficient sites of dynamo activity. 

The theory of magnetic field generation and its subsequent evolution in young PNSs is primarily put forth by TD93 and then 
further elaborated by TMOl . A summary of their findings, relevant for this work, is given below. 

2.1.1. Case of Rapid Rotation-Magnetars 

Based on theoretical and phenomenological arguments, Duncan and Thompson (1992) and TD93 conclude that PNSs with initial 
roation periods near break up (Pjot ^ 1 ms) efficiently convert the nascent star's rotational shear energy into toroidal and poloidal 
flux via a-Vl dynamo action at the base of the convection zone where the overturn time is similarly ~ 1 ms. They argue that the 
free energy available from differential rotation Eq_ for rapidly rotating (P,ot ~ 1 ms) PNSs is roughly Eq_ ~ 10^~ergs. This enormous 
reservoir of energy could in principle, generate ordered fields as strong as Z? ~ 3 x 10'^ G throughout the star. This led them to 
conclude that PNSs born rotating near break up form dipole fields during the convective Kelvin-Helmholtz phase that are much 
stronger than those of pulsars. These highly magnetized neutron stars with dipole field strengths Z?dip ~ 10'"*- 10'^ G are referred to 
as "magnetars." 

2.1.2. Slow Rotation - Analogy with the Solar Convection Zone 

In the absence of rapid rotation or a large reservoir of helicity, turbulent convection alone may lead to sizable magnetic stresses. 
There is some phenomenological evidence for this idea, most of which comes from the Sun. Between 1986 and 1988, a period 
of increasing solar magnetic activity, the solar 5-minute acoustic oscillations displayed a significant (one part in 10'*) increase in 
frequency (Libbrecht & Woodard 1990). A rise in the rms magnetic stress to values of B^^s ~ 200G at the top of the convection 
zone, and extending downwards to substantial depth with greater strength, could explain the observed trend in acoustic frequency 
(Goldreich et. al 1991). At the top of the solar convection zone, the overturn time (^5-10 min) is much shorter than the rotation 
period (~ 22 days). Therefore, it seems that convective motions alone may sustain near equipartition (with the turbulence) rms fields 
on timescales that are many orders of magnitude in excess of the convective overturn time. 

In PNSs, unlike the solar convection zone, the gradient of the turbulent diffusivity lies in the opposite direction with respect to the 
buoyancy force acting on a given flux rope in thermal equilibrium with its surrounding (Tobias et. al 1998, 2001 ; TMOl). As a result, 
both buoyancy and turbulent diffusion tend to transport magnetic structures of scale Ip, the gas pressure scale height, to the top of the 
convection zone approximately in an overturn time.^ At the same time, flux ropes that are generated and amplified by the stretching 
action of the convective cells are pinned at the base of the overshoot region. 

One of the most commonly held beliefs in dynamo theory is that turbulent magnetic fields generated and amplified by convective 
overturn come roughly into equipartition with the turbulent kinetic energy. Mixing length theory provides for a rough estimate of the 
convective velocities at the base of the convection zone 

F, ^pv]^ V, ~ 2 X \(fFlJ^ p['f cms-' (2) 

where Fc and Vc are the convective energy flux and velocity amplitude, respectively. This yields a saturation field strength 

fisat ~ V4^vc ~ 6 X 10'^ plf G (3) 

on the scale of a pressure scale height Ip at the base of the convection zone where Ip ^ I- 10 km. If a near-equipartition magnetic flux 
tube generated at the the base of the convection zone is "pumped" and pinned to the top of the PNS convection zone in the overshoot 
region (KJM96, TMOl), then the local field strength may approach equipartition with the matter near the neutrinosphere for typical 
values of temperature and density (T^ ~ a few MeV and < 10'~gcm"^) at relatively small depths. 

2.1.3. Long-Lived Large Scale Magnetic Structures in the Limit of Slow Rotation 

The ratio of the length scale of the energy bearing eddies to the radius of the neutrinosphere, Ip/Rv ~ 10~\ is ~ 10^ times larger 
than the analogous ratio for the Sun. This leads one to believe that the magnetic structures that are produced by convective turbulence 
in PNSs are much larger than those of the Sun as a fraction of the total stellar radius (TD93). The value of the fluctuating rms magnetic 
field is given by 

^i-ms ~ ey^Bsat, (4) 

where is the convective efficiency parameter In the Sun, eg ~ 0. 1 (Goldreich et al. 1991) and if we naively extrapolate this value to 
PNSs, we obtain with the help of eq. (3) and with all things being equal, B„as ~ 2 x 10'^. One suspects that these magnetic structures, 

* This is true so long as magnetic flux is passively mixed by the turbulence. 
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with spatial dimension roughly equal to the depth of the radiative layer above the convection zone, are created and destroyed on the 
eddy turnover time ^ 1ms. However, Cattaneo and Vainshtein (1991) point out that larger scale structures with smaller flux densities 
can effectively resist eddy diffusivity for a significant number of overturns, implying a breakdown of kinetic approximation. If a 
starspot possesses a dimension /spot ~ 3/p and is due to the random orientation of the smaller dipoles created on the scale of the 
energy bearing eddies, its field strength Bspot ~ 1 /lOBims ~ 2(6^/0. 1)'^^ x lO'^^G. Such a starspot covers ~ a few % of the surface 
of the PNS and may exist for several hundred turnover times before eddy diffusivity can efficiently transfer flux from large scale 
structures to the eddy scale (Cattaneo & Vainshtein 1991). 

From an even more phenomenological perspective, supergranulation on the surface of the Sun may provide a clue to the mechanism 
of large scale field generation in slowly rotating PNSs. In the Sun's quiet photospheric network, magnetic flux is generated and 
destroyed by granule and supergranule motions, rather than originating from the strongly magnetized active regions (Schrijver et al. 
1997). The intermittent magnetic fields associated with supergranule motion typically approach equipartition with the supergranule 
flow while lasting only a few supergranule lifetimes. 

The characteristic scale of a solar supergranule is ~ 30 times larger than the gas pressure scale height at the photosphere or 
equivalently, ~ 30 times larger than a characteristic granule (Stenflo 1989). Also, supergranules overturn with velocities that are ~ a 
few times slower than the velocities associated with granular motion, which are ~ 1 kms"' in the Sun. If we savagely map the solar 
values for granule and supergranule motion to the surfaces of PNSs, one expects supergranule lifetimes on the surface of the PNSs to 
be anywhere from a few hundreds of milliseconds to a few seconds. The scale height of the neutrinosphere ~ 1/30 the stellar radius 
measured at the neutrinosphere R^, implying that only a single supergranule may encompass a sizable fraction of a PNS. 

So far, we have provided the impetus for studying neutrino bubble instabilities in PNSs by motivating the existence of a small 
number of strongly magnetized starspots. Now, we consider in detail the local dynamics and thermodynamics of their linear driving 
mechanism. The next few sections are quite technical. They are riddled with large algebraic expressions whose meanings may not be 
immediately transparent to the untrained eye. For readers not concerned with the details of linear radiation magnetohydrodynamics 
in degenerate media, we suggest skipping ahead to §5.6 where estimates of stability and driving are given in terms of physical 
parameters appropriate for PNSs. 

3. FUNDAMENTAL ASSUMPTIONS AND EQUILIBRIUM 

We consider fluid perturbations upon a highly conducting magnetized fluid consisting of neutrons, protons, electrons, and photons 
interacting with radiation consisting of neutrinos (and anti-neutrinos) of all flavors. The interaction between normal matter and the 
neutrinos is mediated by the various sources of opacity. We make the simplifying assumption that all neutrino opacities increase with 
the square of the neutrino energy E. In particular, for electron neutrinos and anti-neutrinos, 
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whereas for all other types of neutrinos. 
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(We adopt the same constant fiducial energy scale Eo in both cases.) The prefactors kq and kqx are assumed to be constant and 
independent of composition. Only the electron-type neutrinos can be absorbed and emitted by pair capture processes 

e~ + p ^ n + v^ 

+ p + Ve (7) 

implying that only the electron-type neutrinos can exchange energy with the matter while all neutrinos can exchange momentum 
with the matter via elastic neutrino-nucleon scattering 

^+N^iy+N (8) 

where v represents any of the six (ve, Ve,vx,vx) neutrino species and denotes either neutrons or protons. 
For now, we shall assume that thermal equilibrium 

Tg = T,, (9) 
where Tg is the gas temperature of the gas, and chemical equilibrium 

y.p+Pe = ^'n+t^u.j (10) 

where /i, is the chemical potential of the ;* species, hold in both the equilibrium and the perturbed state^. This assertion is more 
subtle than it appears at first glance and we discuss its accuracy later on in this work. 

Assuming that the neutrinos are optically thick and that the neutron- proton-electron-photons (npe-f) fluid is an excellent conductor, 
the expressions for mass, momentum, lepton fraction, energy, and magnetic flux conservation are 

|^ + V-(pv) = 0, (11) 

^ At temperature and densities typical for young PNSs, strong and electromagnetic interactions easily maintain thermal balance between baryons, electrons, and 
photons. Also since these species are so well coupled to one another, their mean free paths are vanishingly small. For these reasons, from here on we shall refer to 
"the fluid" or "the gas" as the multi-species gas consisting of neutrons, protons, electrons, and photons. 
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and 



V-B = (16) 

respectively. Here, p, v, P, g, B, U , F, Yl, and Fl are the fluid density, fluid momentum, total (n/5e7+neutrinos) pressure, local tidal 
gravitational field, the magnetic field, total energy density, neutrino flux, lepton fraction, and lepton flux, respectively. 

Since PNSs are optically thick to neutrinos of all types, we assume that in eqs. (12), (13), and (14) the diffusion approximation 
accurately describes energy and lepton transport as well as momentum exchange. Thus, we only consider fluid perturbations that 
are optically thick to neutrino diffusion, restricting our analysis to the region below the surface of last scattering i.e., below the 
neutrinosphere. The radiative energy flux F is given by the sum 



F = F^„+F^„+F^^ 

where F, is the energy flux of the /^species. The energy flux for electron-type and X-type neutrinos is given by 
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respectively. Since X-type neutrinos are created in the neutral current pair annihilation process, an excess of /i- or r-type neutrinos 
is absent. As a result, only electron-type neutrinos carry lepton number due to their ability of being absorbed and emitted during pair 
capture and emission processes. The lepton flux is given by 



F, =- 



(20) 



6Tr~KopTr'c~ 

The expressions for neutrino energy and lepton fluxes as well as other important thermodynamic quantities obtained from Fermi- 
Dirac distributions are briefly derived in Appendix A. The degeneracy parameter 77, = fii/kBT will be interchanged with /i, throughout 
this work depending on the particular application involved. The mass density is given in terms of the nucleon or baryon number 
density 

P = «'««, (21) 

where »z„ is the mass of the neutron, implying that the mass difference between the neutron and proton has been ignored when 
evaluating the density. The total pressure and energy density are given by 
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The electron and electron neutrino fractions are given by 
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where the total lepton fraction is defined as Yi = Y^+Y^^. The baryon energy density, Ub, is given by that of a monatomic ideal gas 
near the neutrinosphere, where the effects of degeneracy pressure are negligible for the baryons. However, in the core the zero-point 
Fermi pressure is substantial for the baryons, which is why we have written the baryon equation of state in a somewhat general 
manner since its form depends on the location below the neutrinosphere. 
The condition for chemical equilibrium may be written as 



im„-mp)c^ 



-In 



3/2 



and the fact that the total nucleon number density « is conserved implies 
Finally, charge neutrality implies 



Yp + Y„=\. 



1 + ^ 



(30) 



(31) 



(32) 



This gives us 21 equations for 21 variables; p, v, P, B, [/, F, «, i^, T, 77^, 77^^, Yp, K,,, (F^^ + Fij,), ^vx^ ^L^Jb, (Ue- + Ue+), U-y, 
(U^^ + U^J, andf/^^. 

We treat the background state of the PNS envelope in a fashion similar to that of a normal main sequence star The gravitational 
field is assumed to be fixed, which is a good approximation for our oscillation wavelengths of interest near the neutrinosphere, 
where the radiative flux is roughly constant with depth as well. That is, we adopt a plane parallel approximation. Rather than photons 
transporting energy from one depth to the next, neutrinos of all flavors take on that role. Also, we implicitly assume that chemical and 
thermal equilibrium holds for our static background and perturbations. For PNSs, this condition is met for any dynamical timescale 
of interest below the neutrinosphere (see Appendix C for a detailed discussion regarding chemical and thermal equilibrium). 

The major qualitative difference between the environments of PNSs and the stellar envelopes considered by Blaes & Socrates 
(2003) is that the radiation considered here is degenerate with a large chemical potential. Since the flow is optically thick, gradients 
in both neutrino temperature and chemical potential transport both neutrino energy (not just heat) and particle number For the black 
body (/i = 0) radiation fields considered by Blaes & Socrates (2003), the background radiation flux interacts with short-wavelength 
compressible fluid perturbations in two ways. Radiation piles up downstream density maxima with respect to the background flux 
F, an effect known as "shadowing" (Lucy & White 1980, Lucy 1982). Thus, the radiation flux couples to the density perturbation, 
producing a pressure response cx F • ^Sp. The other effect arises from changes in a fluid element's position during oscillation 
with respect to surfaces of constant background radiation pressure, which results in a pressure response cx ^ • F where ^ is the fluid 
displacement. In the case of PNSs, both F and F^ interact with compressible short wave-length fluid oscillations in this manner. 
Therefore, neutrino bubble instabilities are driven by both F and F£ rather than just F alone. More fundamentally, radiative diffusion 
due to gradients in both temperature and chemical potential are responsible for radiative driving, but otherwise the mechanics of 
photon bubble-like instability is the same as the neutrino bubble instability. 



4. TOTAL PRESSURE PERTURBATION IN THE LIMIT OF RAPID NEUTRINO DIFFUSION 

We wish to study standard magnetosonic motion altered by the presence of radiative diffusion. In the optically thick limit, radiation 
alters the fluid's motion by gradients in its pressure. The instabilities classified by Blaes & Socrates (2003) occur in the hmit where 
radiation diffuses rapidly compared to the local magnetoacoustic crossing time. We consider the same limit here. 

Generally speaking, the total pressure may be written as 

p = p(p,r,^„^,J (33) 

and upon perturbation, we have 
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We may eliminate pe in terms of p, T, and p^, by use of the eqs. (30), (31), and (32) which denote beta equilibrium, baryon 
conservation, and charge conservation, respectively. The total pressure perturbation may then be written as 

As previously mentioned, neutrinos diffuse down gradients of both temperature and neutrino chemical potential. In the limit of 
rapid diffusion, the spatial distribution and relative amplitude of the neutrino temperature and chemical potential are almost purely 
determined by the divergence of the neutrino and lepton flux. To see this clearly, we write the linearized expressions for conservation 
of lepton number, eq. (13), and the first law of thermodynamics, eq. (14). This gives us 



-iujSYi + 5\ ■ VFz, = -iuj 



dYL 
dT 



-iLuSU + S\-'^U + i(U + P)k-S\ = -iu 



+ 5v-WYl = — k-(5FL 
n 



(36) 



+ (5v • Vt/ + /(t/+P)k • iJv = -ik • (5F. 



(37) 
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(5Fl and 5F are given by 



and 



P 



ilC 



5F: 



where 



IC = 



en 



and 



1 



1 



(38) 



(39) 



(40) 



3kox 3ko 

Note that we are operating within the framework of the WKB approximation where all perturbation quantities take on a plane wave- 
like form oc e'C*"-'^') Therefore, our analysis from here on, is only credible for perturbations with wavelengths smaller than the gas 
pressure scale height. 

From eq. (38), the change in lepton flux is determined by gradients in the neutrino chemical potential mediated by the electron 
neutrino opacity as well as a component that arises from changes in density. Eq. (39) implies that changes in the radiative heat 
flux originate from both changes in the temperature gradient mediated by the sum of the electron- and X-type neutrino opacities and 
gradients in the neutrino chemical potential mediated solely by the electron neutrino opacity in addition to a component that comes 
from fluctuations in density. 

A first step in obtaining the total pressure perturbation is to solve for 6p,,y^, in terms of the other thermodynamic variables, 5T and 
5 p. We (arbitrarily) choose the Unearized expression for lepton conservation, eq. (36) as our starting point. We have. 



dT 



5T + 



f dYL 



&P-y,+ 



dp 



6p 



+ «(5v- VyL = ' — k-FL 



It is useful to define a characteristic diffusion frequency, t^L.diff associated with Yl transport 

k^JC 



tt^kqp 



5p^, 



which allows us to write 

{dYjdp)^ 



6pv, = - 



WL,di£f [dYi/dpi,^] 



Sp- 



WL.diff = W + i 



k Fl 



5p 



nujum{dYL/dp.y) ^ p ujudis [dYi/ dpi,) 



-ST- 



S\ ■ VYl 



LOL.dm{dYi/dpi,^ 



(41) 



(42) 



(43) 



BS03 found that in the limit of rapid radiative diffusion, temperature perturbations become reduced on the scale of the mode 
wavelength as long as the characteristic diffusion timescale is much shorter than the mode's period. Here, both the temperature 
and neutrino chemical potential are subject to the same diffusive reduction since neutrino diffusion occurs along gradients of both 
temperature and chemical potential. Following BS03, we recognize that in the rapidly diffusing hmit. 



ipv^ , (5r oc 8p (X k ^5p 



(44) 



where Wdiff represents some characteristic diffusion frequency. Note that we have assumed that u) oc k, which implies that we are 
investigating the properties of short-wavelength MHD waves i.e., the Alfven, fast, and slow waves. Many authors have noted that in 
general, there are two diffusion frequencies associated for both lepton and energy transport in optically thick PNSs (Bruenn & Dineva 
1996, Pons et al. 1999, Miralles et al. 2000, Bruenn, Raley, & Mezzacappa 2004). That is, lepton and energy diffusion together 
evolve down gradients of both lepton fraction and entropy (or alternatively, temperature and neutrino chemical potential). Rather 
than investigating overturn instability phenomena such as Ledoux convection, "neutron-fingers," and semi-convection, we focus on 
questions of acoustic stability in the limit where sound waves (and magnetosonic waves) are highly non-adiabatic. That is, both 
lepton and energy diffusion occur rapidly in comparison to the acoustic crossing time of interest for a given oscillation wavelength. 
Due to our cx E~ opacity parameterization, the "cross-term" for lepton diffusion is absent, constraining lepton evolution to depend 
only on gradients in neutrino chemical potential denoted, here by a single diffusion frequency cjL.diff • 

Now, in the short- wavelength high-fc rapidly diffusing limit, the characteristic frequency associated with lepton diffusion satisfies 

k^JC' 

|wL,diff| ^ ,^ X » c^- (45) 

With this, the neutrino chemical potential is approximately given by 



.TT^Kop 



k^lC 



OYl 
dp 



T.Pi, 



k Fl 

P 



Sp. 



(46) 



The first term on the right hand side produces damping in much the same way as classical Silk damping (Silk, 1968). There are 
however, some subtleties to this concept - a point we will consider later on. The second term cx k • Fl is in part, responsible for 
driving arising from the background lepton flux interacting with changes in the fluid's density. Note that we have made use of the 
fact that \6\\ cx 6p for acoustic and magnetosonic waves, which allowed us to drop the term in eq. (43) cx 6\ ■ VF/,. 
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In order to calculate 5T in terms of 6p, we perform the analysis for rapid diffusion with respect to the first law of thermodynamics. 
We have 




5T + 



dU 



T.P 



= !k F 



dp 
Sp 



U+P 



5p^+5\-[VU-{U+P)V\np] 

2 ^Pv, 



TT^KqP 



(47) 



Now, we may define another characteristic neutrino diffusion frequency, t^E.diff which measures the relative importance of energy 
transport via the diffusion of neutrinos down their temperature gradient 

WE.diff = 1^ + 1 • (48) 

In a similar fashion to 5p^^_, the temperature perturbation becomes 



ST: 



{du/dp,:^ 



WE,diff 



Spur 



{dU/dp)^^^^-{U + P)/p ,dy[VU-iU + P)Vlnp] 



k F 



Sp 



WE.diff (dU/dT)^ 
k^JCp^^, 



■Sp 



WE.diff (dU/dT) 



t^E, 



diff {dU / dT) p TT^KopujE,dis{dU /dT) 



-Spv^. 



(49) 



In the short-wavelength high-fc limit, energy transport by means of neutrino diffusion is rapid compared to the acoustic or magnetoa- 
coustic crossing time such that 

|a;E,d.ffH ^i^. (50) 

Again, by realizing that 5py^_ cx k~^Sp and \5\\ oc Sp for acoustic and magnetosonic waves in short-wavelength limit, the temperature 
perturbation may be written as 



ST: 



. HetfP 
I— 5 — 

k-fCkgT 



/dU_ 



U + P 



^ k-F^ k^lCp,^^ 
Sp-\ Sp+i — T Sp^^ 

P TT^KoP 



and with the help of eq. (46), we have 



ST 



k^lCklT 



dU_ 
dp 



U + P 



- Pu, n 



dp 



k-[F-/i,,FL] 



p 



Sp. 



(51) 



(52) 



T,P., P 

The first three terms inside the braces in eq. (52) represent Silk damping. Upon compression, neutrinos leak out along the temperature 
gradient of the fluid element and cool the gas. Cooling during compression can be thought of as a general criterion for damping. 
The term cx [F-zi^^Fl] is responsible for unstable driving due to a combination of both the energy and lepton flux coupling to 
fluctuations in density. 

In this light, the total pressure perturbation is divided into two parts; one which yields a standard acoustic response oc Sp while the 
other component (x ik~^Sp leads to radiative driving and damping. That is, 

SP=SP„, + SP 



where 



The quantity 



SP,c = -TT- 



dP 

dp 



Sp and 



T,Pu, 



fdP_ 

\'d~p 



T.p.„ 



BP 

Tp 



T.Pu 



dpe 

dp 



T.p„ 



ST + 
dP 

dpe 



( BP 



pJ-V-v 



Spv 



pj 



(53) 
(54) 

(55) 



is the square of the isothermal and "isomeric" sound speed c,-. When pair capture and emission occur prodigiously, the rapid ex- 
change of electron-type neutrinos fixes the gas temperature to that of the neutrinos and the composition of the npe^ gas is determined 
by the evolution of the electron-type neutrino's chemical potential. Due to rapid neutrino diffusion, which depends on the gradients 
of both the neutrino temperature and chemical potential, perturbations in temperature and neutrino chemical potential are smaller by 
a factor cx in the short-wavelength limit. Upon compression, both the temperature and relative chemical composition of a given 
fluid element remain fixed to lowest 0(^"') while its density and electron chemical potential evolve. With this in mind, by combing 
eqs. (35), (46), and (52) the total pressure perturbation may be written as 



SP = crSp + i( — 



k^KklT 



+ i 



dP 

dpv. 



pT 



dU_ 

dp 

dYL 
dp 



U+P 



T,p^ 



T.p„ 




Pu,n 



dp 



T.Pi, 



k- [F-/vFl] 



' Sp 



(56) 



* By "isomeric" we mean that the composition of a fluid element is fixed. However, this is somewhat misleading. Only when the chemical potentials of all of the 
individual particle species are fixed can the composition be thought of as being frozen in. In the limit of rapid electron-type neutrino diffusion, the chemical potential 
of the neutrino is fixed over a wavelength to lowest 0(k^^). Therefore, the relative concentrations of electrons, protons, and neutrons may evolve subject to the 
constraint S/ip + Sfie-SfJ-n = (5/i^,=0. 
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4. 1 . A Note on The First Law of Thermodynamics in the Limit of Rapid Neutrino Diffusion 
Here, we resort to the first law of thermodynamics 

dE = TdS - PdV + Hi dNi (57) 

in order to gain insight into the results accumulated thus far It is convenient to write the first law in terms of U = E /V and Yl rather 
than E and A^, i.e., 

U + P 

dU dp = nT ds + njiy dYi, (58) 

P 

where s is the entropy per baryon. By differentiating with respect to density 

nT(—\ =(—\ -^-^Ji,n(^\ (59) 

we see that the right hand side yields the thermodynamic motivation for Silk damping, when compared to the appropriate terms in 
eqs. (52) and (56). If pressure leads density in time such that 5P cx -iSp, the mode will be damped. Thus, the criteria for Silk 
damping in PNSs is given by 

fds\ fdU\ U+P fdVA 

If the entropy of the fluid decreases upon an increase in density, then the fluid cools when compressed and is damped out (for now 
we are only considering the homogeneous case where the effects of the background F and Fl are small). Note however, that 

dU\ U + P fdYA 

* <0 and — ^ <0. (61) 



^P / 7-,A.„„ P \9P y r,,,,... 

Hence, upon compression it is possible for radiative diffusion to heat the fluid rather than to cool it, making it possible for radiative 
diffusion to drive compressible perturbations unstable. In the fj,^^ — > limit, cooling necessarily accompanies compression as in the 
case of classical Silk damping, where neutrinos only carry heat along with them as they diffuse down a temperature gradient. To 
show this, we again exploit the first law 

dU U+P dp ds dYL 

=-V-t=«r \-np.y_ (62) 

dt p dt dt ' dt 

where F denotes the flux of energy. When the chemical potential of the radiation is zero, the energy flux mediates the transport of 
heat alone. To further illustrate this point, we re-write the first law as 

ds 

«r-=-V- [F-m.FlI-Fl-V^,. (63) 

Burrows, Mazurek, and Lattimer (1981) noted that the second term on the right hand side is always positive for global cooling models 
of PNSs and therefore, always leads to heating as neutrinos escape. Those authors refer to to this effect as "Joule heating" in analogy 
with the well known effect in electrodynamics. ^ Upon perturbation, the first law becomes 

-iujnT[Ss + ^-Vs] = -ik ■ [SF-fi^JFi] +Sfi^y ■ Fl, (64) 

where last term on the right hand side vanishes since Fl is a constant for our equilibrium. This leaves us with 

-iw«rfc~-/k- [S¥-fi^^6Fi,] =-/k-5FQ (65) 

where 6Fq can be thought of as a heat flux and is responsible for cooling the perturbations via radiative diffusion. From eqs. (17)- 
(20), the heat flux Fq = F-pi^¥i^ only depends on gradients in temperature i.e., heat is transported via changes in temperature alone. 
Burrows & Lattimer (1986) point out that this is due to our parameterization of opacity, which depends on the square of energy rather 
than on energy and composition. Thus, the temperature perturbation in the limit of rapid neutrino diffusion given by eq. (52) can be 
written in a more compact, yet descriptive form 



6T~- 



k^ICklT 



LunT \ — ] + 



dp 



Sp. (66) 



Compare eq. (46) with eq. (66). Clearly, the entropy content of a fluid element is tied to its temperature perturbation while its 
lepton fraction is tied to fluctuations in the neutrino chemical potential. This reinforces the concept that heat flows down gradients in 
temperature while degeneracy energy is transported due to gradients in neutrino chemical potential. 

' This point is made clear by realizing that Fl ■ V/ii/, = -|Fl|"/o" where a is the radiative conductivity associated with lepton diffusion and that the lepton flux Fl 
can be thought of as a current density. 
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4.2. The Possibility of Diffusive Heating of Compressible Waves 

In the limit of rapid neutrino diffusion, if the entropy of a compressed fluid element increases appreciably to counter the loss of 
degeneracy pressure that accompanies lepton transport, then the fluctuation gains energy during an oscillation. The criterion for this 
to occur and its applicability to PNSs is sketched out in Appendix B. 

The entropy of a fluid element increases if the baryons exhibit some level of degeneracy as in the case of PNS cores. In order to 
illustrate this property, compare the "Silk determinant" for an ideal non-relativistic, non-degenerate gas of baryons {P cx pT) with 
that of a completely degenerate non-interacting and non-relativistic gas of baryons 

a^ideal dU U+P P a^deg „ 

— — = = — < while = 0. (67) 

op op p p Op 

The fact that the entropy response of a degenerate gas vanishes follows directly from the third law of thermodynamics. In other 
words, a purely degenerate gas is adiabatic hy fiat since it contains zero entropy. Furthermore, in a mixture of baryons consisting of 
protons and neutrons in the core of PNSs, the total equation of state may be super-adiabatic in the sense that the net entropy increases 
upon compression. Consider the zero-point energy density of the neutrons given by eq. (B3). The estimate given in Appendix B 
indicates that upon compression, the entropy of the neutrons increases in the limit of rapid neutrino diffusion. An increase in density 
necessarily leads to a decrease in lepton fraction Yi due to overall charge neutrality. The removal of a proton, which results from 
deleptonization, and the addition of a neutron results in a net increase in pressure or energy per unit volume upon compression since 
the chemical potential i.e., Fermi energy, of the neutrons is greater than that of the protons. 

Another positive contribution to the Silk determinant originates from changes in electron pressure or energy density resulting from 
changes in density. Though the electron pressure does not explicitly depend on fluid density, it implicitly does through the electron 
chemical potential and the condition for chemical equilibrium. In the limit of large electron degeneracy, the electron pressure 
Ue = 3Pe oc p^. Therefore, changes in Ue result from changes in /ig. The condition for chemical equilibrium, eq. (30), informs us that 
Pe increases by a factor of kBT{Yp/Y„) upon compression. Eq. (B6) displays this point. 

The deleptonization of a fluid element contributes to the Silk determinant, quantified by the term cx p^^^ dYi j dp at fixed temperature 
T and neutrino chemical potential p^^ . It is tempting to equate the phenomenon of heating due to deleptonization for compressible 
waves to Joule heating. Such a comparison would be misleading at best. Rather than an energy deposition mechanism cx |F/^|-, this is 
a purely linear phenomenon. The Silk determinant, which is cx [ds j dp) at constant T and v^, denotes only the entropy evolution of a 
fluid element. The deleptonization term merely acts to enforce that the heat content of the neutrinos, rather than heat and degeneracy 
energy, is being taken into account as the entropy of a fluid element changes during compression. 

The sum of all these actions together must overcome the loss of degeneracy pressure accompanied by deleptonization. Damping 
due to loss of lepton pressure is provided by the component of the pressure perturbation that originates from changes in chemical 
potential. This effect occurs in global cooling models of PNSs as well (Burrows, Mazurek, & Lattimer 1981). Not only is thermal 
energy altered by neutrino diffusion, but the large lepton degeneracy energy (particularly in the core) is lost as well. Loss of lepton 
pressure resulting from neutrino diffusion always damps out compressible motion. 

5. COUPLING OF MAGNETOSONIC MOTION TO THE BACKGROUND NEUTRINO FLUX: NEUTRINO BUBBLES 

Below we derive the stability criteria and growth rates for the fast and slow neutrino bubble instability while explaining the nature 
of the driving mechanism as well. Two distinct, yet equivalent, interpretations of these radiative magnetoacoustic instabilities will 
be presented. One comes from understanding the relationship between the radiation pressure and density perturbations. If radiation 
pressure lags behind density (in time), heating accompanies compression - the necessary criterion for driving. In the frame of the fluid 
element, the radiation pressure perturbation is composed of three distinct components. The balance between the three components 
ultimately determines the stability properties of the magnetoacoustic mode. Another way of understanding radiative magnetoacoustic 
driving is by examining the relationship between the fluid velocity (or displacement) and the perturbed radiation force. The radiation 
force perpendicular to the wave vector is given by the sum of the perpendicular components of the linear energy and lepton flux 
perturbation. This drives magnetosonic motion unstable since the fluid oscillates in part along the magnetic field line. That is, the 
uniform background magnetic field B may in general, lie in any direction with respect to the wave vector k, allowing the motion of 
the fluid and the radiative driving force to coincide, leading to a net over-stability. 

We first identify the basic properties of standard magnetosonic waves in the absence of stratification and radiative diffusion. In 
order to understand neutrino bubble instability from the relationship between pressure and density perturbations, we calculate the 
growth rates by use of the magnetosonic wave equation and isolate the relevant components of the pressure perturbation that lead to 
damping and driving. We then perform the exact same analysis on the linearized equation of motion where the component of the 
perturbed radiative energy and lepton flux perpendicular to the wave vector drives magnetosonic modes unstable. 

5.1. Basic Magnetosonic Waves 
Linearizing the continuity, momentum, and induction equations, as well as Gauss' Law, yields 

-iujSp+ipk - S\ + S\ ■'Vp = 0, (68) 

-iujpSy = -ikSP+gSp+ — (k X (5B) X B, (69) 

-iuj6B = ik X (5\ X B), (70) 

and 

k-(5B = 0, (71) 
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respectively. 

The procedure for obtaining the asymptotic growth rate for radiatively driven magnetoacoustic instabilities was carried out in 
detail by BS03. Based on analysis of the general (and cumbersome) dispersion relation, they concluded that in the short wave-length 
limit, the radiatively driven MHD instabilities behaved as the standard fast and slow magnetoacoustic modes to lowest order where 
the effects of radiation produced a relatively small imaginary correction to the mode frequency. In order to isolate the mechanism for 
over-stability, BS03 then performed an asymptotic modal analysis where they solved the momentum equation by using the standard 
fast and slow eigenvectors and applied them to all of the 0{k~^) terms in the momentum equation; the pressure driving terms cx ST, 
buoyancy force cx gSp, and to the Eulerian component of the the density perturbation which is oc (5v • Vlnp. This technique is 
mathematically similar to the work integral approach in classic stellar pulsation theory (Unno 1989) for obtaining growth rates due 
to non-adiabatic driving of nearly adiabatic pulsations. 

To lowest order, the linear momentum equation is given by 



-ibjp5\ = -;k 



B (5B 

47r 



+ 1- 



k B 

47r 



(5B 



where (5Pac = cl6p is the short-wavelength acoustic pressure response to a density perturbation 5p. Also to lowest order. 
Eliminating JB, we obtain 



5p = p ^ and i5B= — (k-(5v)-^-^Jv. 

LO UJ UJ 



-ip — ov = -ik 



[ci + v-f\bp-ci 



2(k-VA)' 



bp 



+ !Xk- va)va bp. 



(72) 



(73) 



(74) 



where va = B/ \J\'Kp is the Alfven velocity and o)^ = - (k • Va)^ . In the optically thick limit, neutrinos force the gas via a gradient in 
their pressure which is || to k in the short wave-length limit. This requires the waves to possess some longitudinal, or compressible, 
component which automatically rules out the purely incompressible shear Alfven wave as a candidate for over-stable radiative 
driving. In order to isolate the magnetosonic fast and slow waves, let's define a mode polarization e(w, k") such that 8y = e{oj, k)tp{uj, k) 
where is some complex-valued amplitude. Choosing ip = Sp/p gives us the magnetosonic eigenvectors. 



— rcf + vi 



k-(k- va)va 



where the two magnetosonic modes must satisfy the dispersion relation 

,~,2 



:k^cj 



^ ,1 ' 

^2 



(75) 



(76) 



In principle, the fast and slow modes posses components both || and _L to k. This mixed longitudinal and transverse behavior is 
central to the MHD radiative driving mechanism (BS03), a point we will now discuss. 



5.2. The Magnetosonic Wave Equation Subject to Rapid Radiative Diffusion: Asymptotic Growth Rate 

Taking the divergence of eq. (69), immediately allows us to restrict our analysis to compressible short-wavelength motions i.e., 
the modes that can be driven unstable by the effects of radiative diffusion. We have 



-iojk ■ S\ = -ik^ 



SP B SB 

p A-Kp 



+ (k.g)^ 

P 



where we have made use of eq. (71). Eliminating i5B using eq. (70) gives us 

\2 



-iuj (k • 5\) = -ik 



SP Ak-Sy) (k-VA)"(5P .(g-VA)(k-VA)^p 
— + vl , I 

p UJ UJ~ p LO^ p 



+ (k.g)^. 

P 



ill) 



(78) 



Dividing the pressure perturbation into its acoustic (5Pac and driving/damping SP contributions and by utiUzing eq. (68), the magne- 
tosonic wave equation becomes 



~,2 



; 2 2 

-k 



S_p 
P 



-iujS\ ■ V\Tip + k 



UJ" SP ^ Sy 

—5- — + iv^ Vhip- 

UJ^ p UJ 



.(g-VA)(k-VA) Sp 1 Sp 
I S+i(k-g)— . 



(79) 



Neglecting terms depending on background gradients or those originating from radiative diffusion leaves us with only the left hand 
side of eq. (79) i.e., the magnetosonic wave equation for a uniform background. The terms on the right hand side of eq. (79) produce 
C(^~') corrections to a magnetosonic wave and are responsible for the growth and damping rates which can be viewed as C(^~') 
corrections to the mode frequency. In order to evaluate the frequency correction arising from radiative driving, we let uj — > w+Aq, 
where Aq is the damping or driving rate which is independent of the wavenumber k. 
It is useful to recast the gravitational acceleration, g, by using hydrostatic balance 



g = 



'-VP- 

P 



1 

c;-Vlnp+- 

P 



dP 
df 



■7^ VT + 



dP 



(80) 



13 



This, along with the magnetosonic polarization vector given by eq. (75) allows us to eliminate all of the terms cx Vlnp in eq. (79). 
The asymptotic prescription for the eigenvalues mentioned above allows us to convert eq. (79) into an expression for the damping or 
driving rate. After a bit of algebra, we have 



uj P 1 P 

dP 
df 



Cj^ 5P .(k-VA) 



+ i- 



-va 

uj-p 

VT + 



dT 
/ dP 



( dP 



Sp 
P 



5p 
P ' 



If we expand 6P by use of eq. (56) while noting that 



vr = - 



ICklT 



[F-^^^Fl] and Vp^^ 



Fl, 



(81) 



(82) 



the growth rate Aq from eq. (81) is given by 



TT^Ko / dP 



KklT V dT 



dP 



pu?' 



dU_ 

dp 



P-T 

U + P 



dYL 

dp 



+ ^]i:I^(kxvA)-(kxFL) 



-Pi^n 



T.pi, 



p.p 

dp 



T.p^ 



+ ^]LZ^(kxvA)-(kx [F-m.Fl]) 



(83) 



5.3. Analogy with the Work Integral 

BS03 noticed that when the components of the radiation pressure perturbation originating from gradients in background quantities 
are non-zero, radiative driving may then occur. Consider the Lagrangian perturbation of the total neutrino pressure defined as 

AP^. = A + + + P^J = I . VP^ + 6P^, (84) 

where ^ = iSv/uj is the fluid's Lagrangian displacement. More explicitly. 



AP„ = 



dP, 
dT 



dP, 



fdp^ 

\dT 



5T + 



f dP, 
\dpu. 



Spy, 



(85) 



The density dependence in the pressure derivatives is discarded since Py depends only on temperature and neutrino chemical potential. 
By using eqs. (46), (52), and (82), the driving component of APj, becomes 



AP^ 



dPy 
dT 

dP^ 
dT 



P<^c 



HeffP 

KklTk^ 



JCklTk 



k^^-[F- /I, Fl] - k • [F - Fl] (k • 

{ dPy 



dPy 



1T~Ki)P 



3<^(kxe-(kx [F-m.Fl]) 



^ dpu, ) J 
ttkqp 



fc2^.FL-k.FL(k.e 



\dpv,)j ICk^ 



(kxe-(kxFL)y, 



(86) 



and by using the magnetosonic eigenvector given by eq. (75), the driving component becomes 



^pdrive ^ • 



.(k-VA) fdP, 
dT 



JCklTk^ 



(kxvA)-(kx [F-a^^Fl]) 



P 



+ 1 



(k ■ Va) / dP„ \ tt-kqp 



(kxvA)-(kxFL) 



P ' 



(87) 



If the radiation pressure maximum is attained after the fluid's density maximum, there will be a net amount of work during the 
compression cycle. Thus, the criterion for instability via radiative driving can be written as 

Im[Ap*AP^] >0 (88) 

where Ap = 5p+£, - Vp is the Lagrangian density perturbation. It can be checked that upon substitution of eq. (87) into the above 
expression, the same criteria over-stability is obtained as that given by eq. (83). 

To gain a deeper understanding of the driving mechanism, consider purely hydrodynamic acoustic motion where B = va = such 
that £ II k. Substituting an appropriate sound wave eigenvector into eq. (86), we see that 



AP'!"^'= cx 



(dK 
\ dT 



KeffP 



P<'e 



JCklTk^ 



^U'k- [F-M. Fl] -k- [F-^, Fl] (k-k) U 



■K~HOP 



dPu 
dpv, J J ICk^ 



;t"k-FL-k-FL(k-k) } =0.(89) 



That is, for purely hydrodynamic sound waves, radiative driving vanishes due to an exact cancellation of two unique effects, which 
we describe below. 

Regardless of the direction in which a fluid oscillates, the background radiative fluxes F and F^ always force a fluid element if 
there is a change in density. We call this effect "shadowing," following the terminology of the O star line-driven wind literature (Lucy 
& White 1980, Lucy 1982). For the sake of exposition, let's ignore lepton number and assume that Fl = Pu, = 0. The component 
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of the Lagrangian pressure perturbation due to shadowing originates from the V • SF term in the first law of thermodynamics and is 
present in eq. (56) i.e., 

'dP^\ Keff .[ dPy\'k-VT 5p __(dP^\'k-VT 



^ pshadow _ 



-(k-a. 



(90) 



-B^ / K- P 

Radiation piles up before the density maximum as the local decrease in mean free path temporarily hinders the transmission of quanta 
(Baker & Kippenhahn 1962). However, the total extinction over a wavelength remains the same since a fluid element's mass is a 
constant during the compression cycle. In order to conserve quanta, a deficit of radiation, or a shadow, appears upstream of the 
density maximum. The surplus of radiation downstream and the deficit of radiation upstream of the density maximum induces a 
pressure differential that peaks at the point of maximum density along the wave. 

The contribution to the Lagrangian radiation pressure perturbation due to changes in the displacement vector ^ is given by 



^pdisplE 



(91) 



When a fluid element moves upward along the direction of F, the background pressure it samples is smaller than the value at its 
original position. Conversely, when a fluid element moves downward, the radiation pressure perturbation increases. As previously 
mentioned, if ^ j| k as in the case of a purely hydrodynamic sound waves, the pressure perturbations from shadowing and motion 
along the temperature gradient cancel one another out, a point clearly seen when comparing eqs. (90) and (91). 

This harmony between AP^''^''°™ and /\p^'^P''"=<= jg broken when the magnetic field is introduced. The change in radiation pressure 
due to a change in position does not in general match the acoustic response followed by a change in density since the fluid element 
oscillates in part along the direction of the equilibrium B. 

5.4. Lunit ofBlaes & Socrates 2003 

The instability mechanism considered in BS03 is essentially the same as in this work. The major difference was their assumption 
of zero chemical potential for the photon radiation field as they were considering backgrounds that were close to perfect black bodies. 
Setting = {dP/dfi^^) = in eq. (83) we have 
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which is completely equivalent to the one-temperature MHD growth rates given by eqs. (93) and (107) of BS03. 



5.5. Magnetosonic Waves Forced by Radiative Diffusion: The Perpendicular Radiation Force 

An alternative, and more visual, physical picture for these instabilities arises from considering forces on oscillating fluid elements. 
For this particular application, it is both helpful and revealing to adopt the Lagrangian displacement ^ as our fundamental fluid 
variable. Other important fluid quantities are related to ^ in the following way 

^=-/(k.e-|-Vlnp, and ^ 

P P 

This allows us to write the Euler equation, eq. (69) in terms of Lagrangian displacements 



(5B = V X (^ X B) = /k X (^ X B) , 



-/cf(k.a-c?^-Vlnp. 



(93) 



^~-(k.VA)'|-k(c2 + vi)(k.a + VA(k-VA)(k.0 + k(k.VA)(|-VA)-/k 



dt 

where SP = SPdme + SPd^mp and from eq. (56) 
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The lowest order acoustic and magnetic tension restoring force terms on the right hand side of eq. (94) can be reduced to 

(k • va)' 1+ k (c2 + v2 ) (k . a - va (k • va) (k • I) - k (k • va) (I • va) = k^vl^^ (96) 

by use of eqs. (75) and (76) where Vpj, = /k^ is the square of the phase velocity. The local gravitational acceleration g can be 
expressed in terms of the radiative and lepton fluxes by use of hydrostatic balance given by eq. (80) and eq. (82) such that 

' dP\ Keff r„ „ 1 f dP \ TT-Ko, 
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Inserting this into the equation of motion yields 
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The driving terms on the right hand side of eq. (98) may be recast in the following form 
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where the last term is orthogonal to to fluid displacements and therefore cannot drive the mode. Ultimately, the linear radiation force 
per unit mass STr^d is given by a gradient in the total neutrino pressure = V [f^^ + and upon perturbation 



ST,,d=-VP,^-i-SP,. 
P P P 



(100) 



where Py = P^^ +Pv, +Pvx- Following the discussion surrounding eq. (65), the first term may be broken into a component which 
carries heat and one that carries lepton number 
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P P \oT 
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where we have kept the leading terms 5p, the Unear driving force (5J-drive can be written in terms of the two perturbed fluxes 

'dP^ 
dT 
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Thus, the component of the perturbed heat and lepton flux ± to k leads to secular driving of the magnetoacoustic waves (see Figure 2). 
Due to magnetic tension, the fast and slow magnetoacoustic waves possess polarizations which have components along B, which in 
general, are not \\ to k. Note that (5J-drive is the driving force originating from radiation pressure alone. As previously argued, changes 




Fig. 2. — Geometry of the neutrino bubble instability. Here, (5F^ is any general perpendicular radiative flux, which denotes the direction of both (5Fl and <5Fq 
- assuming that the equiUbrium Fq and Fl lie in the same direction. For the sake of clarity, the polarization properties of the depicted wave are that of the slow 
magnetosonic wave in the limit that v\ S> c^. As a result of this choice, the velocity perturbation is almost perfectly || to the equilibrium magnetic field. The driving 
force is || SF± and possesses a finite positive projection with the velocity perturbation 5v, allowing for the stratified background radiation field to perform work on 
the oscillating fluid. 
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in pressure due to temperature and neutrino chemical potential perturbations contribute to driving. Below the neutrinosphere, weak 
interactions occur so rapidly that the temperature and composition of the gas are dictated by the diffusive behavior of the neutrinos. 
With the help of eq. (99), the total driving force becomes 



5T, 



total — ~ 



(103) 



Both (5Fq and (5Fl are cx Sp/p implying that the perpendicular radiative fluxes reach their maximum possible upward values during 
density minima. Upon expansion radiation from the equilibrium heat and lepton fluxes more readily leaks out perpendicular to the 
direction of wavefronts as to avoid regions of enhanced density. Due to magnetic tension, then fluid oscillates in part along B, which 
in general has some projection _L to the wavefronts. Thus, the radiative driving force and the fluid velocity may oscillate in phase for 
certain propagation directions, ultimately allowing driving to occur. 

5.6. Estimates of Local Stability and Driving 

We have shown that in the limit of rapid neutrino diffusion, magnetoacoustic waves can be radiatively driven by both a background 
Fq and Fl. Likewise, the waves can be damped by radiative diffusion due to the conversion of work into heat when {ds/dp)T.,j, < 
as well as the loss of lepton pressure as a parcel of gas deleptonizes, which occurs as long as {dYi/ dp)T.^,, < 0. Clearly, the growth 
rate given by eq. (83) can be thought of as two distinct stability criteria, one arising from the transport of heat and one due to the 
transport of lepton number Ignoring all geometric factors (which are not completely negligible) the "condition" for radiative driving 
due to heat transfer is roughly given by 

^e>^n^4nc„ (104) 

while the condition for radiative driving which derives from lepton transport reads 

- 2 3 

FL>^nYLC^^nYLCi. (105) 

The above stability criteria apply for the slow mode (which always has a larger growth rate than the fast mode; BS03) in the limit 
that both /cj < 1 and the background pressure support is primarily provided by the entropy of the baryons, which happens to be 
the case near the neutrinosphere. The two conditions for instability can be reconciled by re-writing eq. (105) as 

3 3 

kBTFL~YLFQ>^ n kB T Y^ a ^ FQ>'-^PhCi = pv\. ( 1 06) 

Therefore, neutrino-driven magnetoacoustic instability can be summed up by this criterion for both transport processes. 
Substituting input conditions appropriate for regions near the neutrinosphere yields the instability condition 

INSTABILITY ( \ f ^ V ( \ f > 10"' (107) 

VlO^ergs-iy \ Ry J V10"gcm-V \ B J ^ 

where Li,, R^, andp^ are the neutrino luminosity, radius of the neutrinosphere, and density at the neutrinosphere, respectively. The 

stability criteria given by eq. (83) indicates that for field strengths at and above the equipartition value such that B^/87r ~ Ph, the 

Alfven speed in the above stability criteria can be replaced by the isothermal gas sound speed c, ^ (Ph/ pYl^. The equipartition value 

of the magnetic field Beq near the neutrinosphere is roughly given by 

/ ^ 1/2 / y, X 1/2 
Beq =i 2 X 10'^ I — -r^ A { ) G, (108) 

^ VlOiigcm-V \AMeN J 

which can be thought of as the magnetic field's limiting value for eq. (107) to be valid. Deeper in the star, the degeneracy pressure of 
both the baryons and electrons overwhelms the contribution from the baryons' entropy and the Fermi energies of both the electrons 
and baryons, rather than temperature, will enter into eq. (108) yielding a larger sound speed and thus a larger Beq- 

Though the stability criteria for both lepton and entropy driving are roughly equivalent, their respective growth and damping rates 
are not. From eq. (83), we see that the driving rate 7/, originating from gradients in p^^, which are cx Fl, for the slow mode is 
approximated by 

t^FL=-—( — \ \Vp.\ 

(109) 

2 a P Ci 

The factor Pl/P, where Pi is the lepton pressure, may not be entirely negligible near the neutrinosphere. The driving rate due to 
gradients in temperature 73, mediated by Fq, can be estimated in a similar fashion 

(110) 

2 Ci P C: 
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where Pb is the baryon pressure. Near the neutrinosphere, Pb ^ P, but P^ is certainly non-negligible. Noting that P^ ~ Pf ^ {rje /4) Yg Pt, 
indicates that if radiative driving is allowed to occur at relatively large depths, where the degenerate lepton pressure becomes 
formidable, then the neutrino bubble instability results primarily from lepton, rather than heat, transport. 

In terms of equilibrium parameters near the neutrinosphere, the combined neutrino bubble growth rate 7nb = 7l + 7q is approxi- 
mately 

,/ B \ /4MeV\ /10"gcm-^^\'^" / g \ , 

As in the case of the instability criteria eq. (107), the value of the magnetic field B in the slow mode growth rate possesses an upper 
limit whose value is given by eq. (108). 

Let us summarize the results of this section. Near the neutrinosphere, the condition for local driving due to the presence of Fq and 
Fl are roughly equivalent. This is true primarily because the neutrino chemical potential /i^^, ~ ksT, which implies that the flux of 
leptons directly maps to the flux of heat modulo a factor that is ~ Yl. Neutrino bubble instabilities can grow on millisecond timescales 
for strong magnetic field strengths, with growth times > 1 ms for equipartition magnetic fields. Driving due to the presence of Fq 
mostly likely dominates the driving near the neutrinosphere over the radiative driving provided by Fl. 



6. MONTE CARLO CALCULATION 

Assuming that radiative driving occurs in the upper layers of a PNS, neutrino bubbles may then lead to substantial luminosity 
enhancements and large density fluctuations in strongly magnetized regions. In order to make quantitative statements regarding the 
non-linear outcome of the neutrino bubble instability, a theory of radiation-driven magnetoacoustic turbulence is required. To our 
knowledge, such a theory unfortunately does not exist. We speculate how the saturation amplitude may change for different values of 
the local magnetic field strength in §7.2 under the assumption that significant non-linearities occur even for gas pressure dominated 
equilibria. 

In order to proceed, we approximate the effects of the neutrino bubble instability by artificially altering the neutrino opacity in a 
radial and angular cut surrounding the neutrinosphere at early times after bounce. By doing so, we capture the cumulative effect that 
we expect neutrino bubbles produce on the surface of a PNS. From this deformed stellar background, the emitted neutrino flux is 
simulated via a Monte Carlo calculation so that the level of asymmetry in the outgoing radiation beam can be estimated. The angular 
profile of the neutrino flux allows us to further estimate the magnitude of the instantaneous kick as well as the gravitational wave 
strain and r-process efficiency. 



6.1. Numerical Inputs and Techniques 

There are three ingredients that go into our numerical model of the neutrino bubble powered mechanism: atmospheric structure, 
parameterization of the neutrino bubble instability, and radiative transfer. We describe each of them below. 



6.1.1. Proto-neutron star structure 

The PNS that we use in our simulation was born from a ISM© progenitor (sl5s7b - Woosley & Weaver 1995). Collapse, bounce, 
and early contraction are evolved up to a time of 130 ms after bounce using a 1-dimensional core-collapse code (Herant et al. 1994; 
Fryer et al. 1999). To this point, we assume that the evolution of the star is spherically symmetric and neutrino bubble oscillations 
have not yet begun. 



6.1.2. Neutrino Bubbles as a modified opacity 

In our stellar model, the onset of neutrino bubble instability introduces deviations from spherical symmetry. As discussed in §2. 1, 
neutrino bubbles are driven in starspots above the convection zone, but below the neutrinosphere. Effectively, one can model the 
global radiation asymmetry due to potential luminosity enhancements in starspots by reducing the opacity in those regions while the 
star's hydrostatic structure remains the same. By doing so, optically thick neutrinos preferentially leak out towards these starspots, 
producing a cumulative asymmetry in the neutrino beam, which emulates the effect of the acoustic instability. If in reality, the 
neutrino bubble instability behaves as the photon bubble instability in the non-linear regime, then large density fluctuations occur 
and the radiation preferentially leaks out of the rarefied regions, locally enhancing the outgoing flux (Begelman 2001, Turner et al. 
2004). 1° 

The angular extent of a starspot is parameterized by a covering fraction C and values for C are taken to be from 0.01-0.1. The 
radial extent for a neutrino bubble-active starspot is parameterized by an inner- and outer-radius 7?.,„ and TZom, respectively. A largely 
phenomenological justification for these values was given in §2.1. Throughout the neutrino bubble-active region, we reduce the 
neutrino opacity by a factor/. This reduction in opacity crudely characterizes the potentially multi-dimensional saturated state of the 
neutrino bubble instability into a single parameter. A Ust of values for C, 7?.,„, TZom , and / that are employed in our calculations are 
given in Table 1 . 

Drawing conclusions regarding neutrino bubble saturation in PNSs from work done on photon bubbles in luminous accretion flows is dangerous. Those authors 
considered radiation pressure dominated media with superstrong magnetic fields. For PNSs, the magnetic field is not superstrong and the radiation field is significantly 
sub-Eddington in neutrinos. 
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6.1.3. Neutrino Transport 

We simulate the transport of neutrinos through our neutrino bubble-deformed PNS 130 ms after bounce via the Maverick Monte 
Carlo transport code (Hungerford et al. 2003). In our calculations, we make several simplifying assumptions. We only model the 
electron neutrinos Ve, while assuming that the asymmetry they experience is representative of all neutrino species. We also assume 
that the asymmetry in the neutrino emission remains the same for the entire phase of neutron star cooling, which will inevitably 
lead to an overestimate of the net asymmetry. With these caveats, we can then calculate the total momentum kick and gravitational 
wave signal based on the total energy emitted in neutrinos for our neutron star (for our model, the total energy emitted is roughly 
5 X 10^^^ erg). 

Neutrino/matter interactions are calculated using a simplified opacity given by (Janka 2001) 

_3a^aoJ4)j,^ (112) 



and 



24 m„ /„.,^2\2 



7^^(F„ + Fp). (113) 



Here, ctq = 1.76 x 10"^"^ cm^ is the characteristic weak interaction cross section and a is the vector coupling constant in the 
vacuum. The absorption opacity is due to electron capture and the scattering opacity Hsc results from elastic scattering off of 
nucleons. Neutrino annihilation was not included as an opacity source since the annihilation cross section becomes vanishingly 
small near the neutrinosphere, when compared to other sources of opacity. This follows from the fact that neutrinos propagate almost 
purely radially. 

Neutrino emission rates from electron capture and electron-positron annihilations are taken from the temperature and density 
dependent functions given by Colgate, Herant & Benz (1993) and the energy distribution for the emitted neutrinos is taken from a 
simplified functional fit to the Fermi-Dirac distribution, as described by Keil, Raffelt & Janka (2003). Emitted neutrino energies were 
sampled from their Fermi sea with a two step (inversion followed by rejection) sampling technique. 

Using appropriate Monte Carlo estimators, we track the spatial distrubition of neutrino energy deposition, luminosity with respect 
to polar angle, and the net neutrino momentum. Neutrino deposition is tallied in a 2D cylindrical grid with 88 radial and 176 vertical 
bins. The Monte Carlo estimator for this quantity is the sum of packet luminosity for all packets absorbed at cylindrical coordinate 
(r,z). The neutrino packets have uniform luminosity weight, so the Monte Carlo noise associated with the calculated deposition is 
proportional to the packet luminosity times ^N(r,z), where N(r, z) is the number of packets absorbed at cylindrical coordinate (r,z)- A 
typical value for the noise is roughly 0.003%, though along the polar axis the noise is significantly higher as the sampled cylindrical 
volume drops to zero at small radii. This purely numerical effect can be seen as dips in the deposition profile along the poles (see 



Table 1 



Monte Carlo Simulations 



c 


/ 




7^,■„ 


Vkick 


< IlTT >e 






(km) 


(km) 


(km/s) 


( 10 kpc) 


0.0 


1 






100 


l.OE-20 


0.1 


2 


90 


50 


1675 


1.3E-19 


0.1 


2 


90 


60 


965 




0.1 


1.33 


90 


50 


700 




0.1 


1.33 


90 


60 


520 




0.1 


1.11 


90 


60 


350 




0.1 


1.11 


90 


60 


220 




0.05 


2 


90 


50 


890 




0.05 


2 


90 


60 


530 


5.E-20 


0.01 


2 


90 


50 


175 




0.01 


2 


90 


60 


120 




0.1 


2 


120 


50 


2350 




0.1 


2 


120 


60 


1500 




0.1 


1.33 


120 


50 


950 




0.1 


1.33 


120 


60 


700 




0.05 


2 


120 


50 


1250 




0.05 


2 


120 


60 


800 




0.01 


2 


120 


50 


250 




0.01 


2 


120 


60 


200 





19 



Energy Deposition [erg/s] 




40 80 120 160 
X [km] 

Fig. 3. — Neutrino energy deposition contours for C = 0.1,/ = 2, TZom = 120 km, 7?.,„ = 50 km. Contour levels plotted are 0.06, 0.2, 0.8, 2.2, 7 and 27 in units 
of lO'^ergss"' . Note that the contours move inward slightly over the cone of lower opacity (cone opening angle is roughly 37° for C = 0.1). This slight deviation 
from spherical symmetry will not likely cause strong variations in the PNS wind. The decrease in deposited energy right along the symmetry axis is the result of 
insufficient sampling due to the smaller volume at small radii. 



Figure 3). Note that numerical noise is given in terms of energy per unit time since we are only considering a short time-slice in the 
evolution of the PNS. 

The angular distribution of the neutrino luminosity is similarly calculated by summing the luminosity of escaping neutrino packets 
into 1 1 angular bins along the polar direction. For a given bin direction 9, packets escaping with direction 6 ± A9 are counted in the 
sum (where A6 = 2.5°). As before, the noise goes as the packet luminosity times ^/N{9), where N{9) is the number of packets in 
angular bin 9, which works out to be roughly 0.01 % (see Figure 3). The Monte Carlo estimator for net neutrino momentum is given, 
in component form {p.^, py and p^), by 




(114) 



where A'fic is the number of escaped packets, Lpch is the weight of the packet in units of luminosity, c is the speed of light and fi^-, fj.y, fj,~ 
are the direction cosines with respect to the coordinate axes. This is the most well determined quantity from these simulations, as it is 
tabulated from all escaping neutrino packets {Ncu- ~ 2 x 10*). The uncertainty in the net neutrino momentum per unit time resulting 
from angular discretization is of order lO^^gcms"^. As the first line of Table 1 shows, numerical noise alone is responsible for 
momentum asymmetries in our calculations, leading to velocity kicks ^ lOOkms"' if this artificial asymmetry where to continuously 
discharge its momentum onto the star for ~ Is. 

6.2. Numerical Results 

Now we are in the position to provide estimates of the kick imparted, r-process yields, and gravitational wave strains which may 
result from the neutrino bubble instability. Values for these observable quantities are provided in Table 1 . 
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Our work focuses on the origin of neutron star kicks, which is the most dramatic effect of neutrino bubble instabihty on the core- 
collapse enviroment. The total momentum carried away by the neutrinos in the collapse is equal to the total neutrino energy emitted 
divided by the speed of light. The asymmetry in neutrino emission carries away a net linear momentum. The total asymmetry is 
measured at infinity rather than at the neutrinosphere. Thus, the optically thin matter above the region of neutrino bubble instability 
reduces the asymmetry by a multiplicative factor ~ e^" and most likely leads to an underestimate of the momentum kick. As 
expected. Table 1 shows that neutrino bubbles may drive kicks in excess of 1000 kms"'. We assume that the kick is discharged for 
~ Is, which can be thought of as the best case scenario. Further, we are implicitly assuming that the spin and the starspot are exactly 
aligned. This assumption overestimates the kick by the value of the cosine of the angle between the spin axis and the starspot. 

In order for neutrino bubbles to affect the r-process or the supernova explosion, it need not alter the total energy deposition of 
neutrinos dramatically, but it must induce an asymmetry in this deposition. Our simulations show that this instability produces only 
a small asymmetry, and the total heating from the proto-neutron star is essentially spherically symmetric, even for our most extreme 
case (see Figure 3). 

However, these small asymmetries may lead to large differences in the gravitational wave output and spectrum. The gravitational 
wave amplitude {hj^) for an observer orthogonal to the symmetry axis for asymmetric neutrino emission is given by (Miiller & Janka 
1997; Fryer etal. 2004): 

2G ^ 

(/jLVe= ^I]Af^(l+z/r)(2xVr2-l)AL,. (115) 

The neutrino luminosity, AL^, at a position x, y, z (r = (x^ +z^y^^) over a given time (At) must be summed over all time to 
give the total wave amplitude. For our suite of models, the gravitational wave amplitude (assuming the asymmetry persists for 
1 s) is sufficiently high to be detectable by advanced LIGO observing a Galactic supernova. ' ' Although this signal will be very 
different than the gravitational wave signal from rotating supernovae, it will not be dramatically different from other asymmetries in 
supernovae (e.g. asymmetric collapse - Burrows & Hayes 1996; Fryer 2004). 

7. UNCERTAINTIES, COMPARISON WITH OTHER KICK MECHANISMS, AND FUTURE WORK 

The neutrino bubble powered kick mechanism presented here is subject to several uncertainties, many of which have not been 
mentioned thus far. In what follows, we assess its major weaknesses and by doing so, provide motivation for future work. 

7.1. Global Stability 

We have focused on the physics of local magnetoacoustic stability in the presence of a stratified degenerate optically thick radiation 
field. We can not stress enough, that drawing conclusions regarding stability from a local linear analysis is dangerous at best. In order 
to confirm that neutrino bubbles are important in PNSs, a global analysis must be performed using a full non-adiabatic treatment 
since the region of greatest driving lies near the neutrinosphere, where the radiative diffusion time is short in comparison to the 
acoustic crossing time. Our kick model requires a magnetic field that perhaps originates from vigorous convective motion beneath 
the upper radiative layer. Turbulent motions almost certainly reduce the driving rate of a given trapped wave (Goldreich & Keeley 
1977). However, even if the interaction between a mode and turbulence is accurately quantified, calculating the cumulative damping 
rate throughout the run of the star requires knowledge of the mode's profile as a function of depth. 

7.2. Saturation Amplitude 

Assuming that one finds trapped neutrino bubble oscillations in a global eigenmode calculation, computing the saturation amplitude 
is a much more difficult task since a fully non-linear degenerate radiation magnetohydrodynamic analysis is required. But, from a 
computational expense point of view, the situation may not be as costly as in the case of main sequence stars and the accretion disks 
that are thought to power active galactic nuclei and X-ray binaries. Since the sound speed is relatively close to the speed of light in 
PNSs, the ratio of the sound crossing time to the radiation diffusion time near the surface of last scattering is only a factor of ^ 30 or 
so. For the A-star discussed in §2, this ratio is ^ 10"*. 

We expect that the saturation amplitude becomes smaller for decreasing values of the magnetic field strength. If the saturation 
amplitude depends on the ratio of equilibrium magnetic to gas pressure evaluated at the neutrinosphere, then large flux enhancements 
would not occur for local field strengths < 10'^ G. However, if the saturation amplitude varies with respect to the ratio of the Alfven 
speed va to the gas sound speed c,, then perhaps field strengths of order ~ a few times 10'"* G may lead to a saturation amplitude 
of order v^/c, ~ 1/3. Of course, the entire kick model relies on the fact that significant luminosity enhancements are reached for 
thermally- and chemically-locked perturbations in gas pressure dominated media. The results of §§4, 5, and of BS03 inform us 
that dynamical growth rates can be achieved for wavelengths of order the gas pressure scale height (where our WKB analysis breaks 
down). Therefore, the energy injected into the mode by the instability is comparable to the kinetic energy of the mode itself. Whether 
or not this leads to order unity saturation amplitudes is a question that can only be answered with confidence by a careful non-linear 
analysis. 

7.3. Magnetic Field Structure 

The neutrino bubble instability can not power a kick unless strong magnetic structures, with length scales significantly larger than 
the gas pressure scale height at the neutrinosphere, are maintained for hundreds of convective overturn times. Reasons for believing 

' ' Depending upon the duration of the explosion, the neutrino luminosity can vary considerably over this 1 s time peiiod and it is possible that the gravitational wave 
amplitude may be even higher (perhaps lower) than this predicted value, but this value gives a rough estimate of the expected signal. 
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that such fields exist in young PNSs were given in §2.1. Our arguments regarding magnetic field structure were primarily based on 
phenomenological arguments derived from the Sun. However, even after the discovery of the solar tachocline (Goode et al. 1991) 
and the understanding of its importance with respect to magnetic field generation (Spiegel & Zahn 1992), the evolution of the Sun's 
surface magnetic field throughout the solar convection zone evades a coherent theoretical picture. Numerical simulations of PNS 
magnetic field structure may prove to be more fruitful than those of their main sequence counterparts since a PNS convection zone 
encompasses relatively few gas pressure scale heights in comparison to say, the solar convection zone. 

7.4. Handling of the Neutrino Transport 

We have ignored the compositional dependence of the neutrino opacities in our linear analysis. Not only does this change the 
numerical value of the opacity, but more importantly, we may have overlooked certain instability phenomena as well. In general, the 
neutrino absorption opacity depends on the proton fraction (see eq. (112) and Janka 2001). Since the proton fraction Yp depends on 
fluid's density, the opacity perturbation 5k varies with the density perturbation 5p implying that standard hydrodynamic sound waves 
which display some degree of radial propagation can be driven unstable by the presence of the radiative fluxes (Unno et al. 1989, 
Glatzel 1994, BS03). Therefore, when the compositional variation of the absorption opacity is taken into account, K-mechanism and 
strange mode phenomena may ensue on the entire surface of the proto neutron star i.e., radiative driving occurs in un-magnetized 
regions as well. If thermal locking holds, the growth rate near the neutrinosphere from local radiative driving is comparable to that 
given by the upper limit of eq. (111). Certainly, if the K-mechanism or strange mode instability occurs at early times, luminosity 
perturbations of only ^ 10% could play an important role with respect to shock ejection and the overall explosion mechanism. 
Fully understanding how this different instability phenomenon affects the neutrino bubble kick mechanism requires answering the 
questions posed by §§7. 1-7.3. 

7.5. Comparison with Other Kick Mechanisms 

If these theoretical hurdles are overcome, a neutrino bubble powered kick mechanism possesses attractive differences when com- 
pared to other mechanisms. In order to obtain spin-kick alignment, as is seen in the Crab and Vela systems, the impulse must be 
imparted on timescales significantly in excess of the rotation period (Lai, Cordes, & Chernoff 2001). In "mass-rocket" models, a 
PNS receives various impulses with durations of order ~ a few milliseconds each (Burrows & Hayes 1996; Lai & Goldreich 2000; 
Fryer 2004) . Therefore, spin-kick alignment most likely does not follow from the hydrodynamic mass-rocket mechanism unless the 
initial rotation period is < a millisecond. However, this fact only rules out the mass-rocket model for the Crab and Vela systems, 
which may not be representative of the Galactic pulsar distribution. Fryer (2004) argues that the mass-rocket works only for rapid 
explosions so that the convective motions in the gain region do not have enough time to damp out the mean motion of the PNS. 
Therefore, mass-rocket kick models are heavily prejudiced towards neutron stars born from progenitors with mass < \2Mq. 

Recently, another hydrodynamic powered kick model that depends on delayed explosions triggering global Rayleigh-Taylor in- 
stability has been further elaborated by way of detailed numerical simulation (Scheck et al. 2004). Their model requires the stellar 
neutrino emission to be somewhat uniformly discharged over the first second so that the refreshment epoch is prolonged. This allows 
global Rayleigh-Taylor type instability to develop beneath the shock front since the density and pressure gradients lie in the opposite 
direction across the shock (Thompson 2000, Herant 1995). A prolonged refreshment is necessary in order to allow the lowest order 
(/ = 1,2) modes to dominate the non-linear phase. This mechanism does not lead to spin-kick alignment and most likely requires 
weak explosions. 

The parity violation model of Arras &. Lai (1999a, 1999b) naturally produces spin-kick alignment as do all neutrino-powered kick 
mechanisms. This results from the basic fact that the Kelvin-Helmholtz time, which is the discharge timescale for neutrino powered 
kicks, is most likely much longer than the rotation period. Also, neutrino-driven kick mechanisms are not subject to convective 
damping in the gain region as in the case of the mass-rocket mechanism since the the shell of accreting material beneath the outgoing 
shock can easily adjust its center of mass to that of the kicked PNS. As previously mentioned, the parity- violation mechanism requires 
global magnetic flux densities ^ 10"' G, a large value even for magnetars. Say for example, that the neutrino bubble powered kick 
mechanism requires flux densities of comparable strength in order to produce order unity luminosity enhancements over a starspot. 
Then, such a starspot would only need to cover ~ 10% of a hemisphere. Thus, the total amount of flux necessary for the neutrino 
bubble mechanism, for this rather extreme case, is ~ 10% of that required for the parity violation model. In terms of the overall 
magnetic energy budget, the neutrino bubble mechanism requires ~ 1% of the total magnetic energy with respect to the parity 
violation model. 

There may not be a single explanation for neutron star kicks. Perhaps all, some, or none of the kick mechanisms mentioned here 
(including ours) work for different classes of progenitors. The best way to differentiate between the various kick mechanisms is to 
constrain models with observations and analysis of the Galactic pulsar and neutron star population. For example, determining the 
fraction of systems which possess spin-kick alignment, quantifying the robustness of the inferred bimodal distribution as well as the 
correlation between kick velocity and magnetic field strength will help constrain the applicability of these various mechanisms. 
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APPENDIX 

APPENDIX A: THERMODYNAMIC QUANTITIES AND RELATIONS 

Here, we derive the expressions for the neutrino energy and lepton flux as well as other important thermodynamic quantities. Much 
of this Appendix is covered by Bludman & Van Riper (1978) as well as Lai & Qian (1998). 
All spin 1/2 particle species are assumed to possess a Fermi-Dirac distribution / where 

1 



fi = 



(Al) 



g(E-tii)lkBT _,_ 1 

for the species. Thus, any given particle distribution is parameterized solely by a temperature and chemical potential, both of 
which vary as a function of depth. In the diffusion approximation the energy density, number density, energy flux, and particle flux 
are given by 



U= / dEUE = 
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E ' 3p, 

respectively. For sake of simplicity, we assume that k{E) depends solely on energy as given by eq. (5) and is the sum of both the 
scattering and absorption opacities. Integrating over energy we obtain 



where F„{r]) is the standard Fermi-Dirac integral given by 
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The chemical potential of the electron-type neutrinos obey the relation p,,^^ = -pu^ which allows us to write 
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and 



Fl = Hjy^ - Hp 



(A8) 



The lepton flux, Fl, is solely carried by the electron- type neutrinos due to their macroscopic mean free paths and because the X-type 
neutrinos are produced only from pair annihilation processes. The above relations were obtained with the help of the following 
expressions (Bludman & Van Riper 1978) 

Fq(v)-Fo(-v) = V 

2 2 
7? TT 

Fi{v)+Fi{-v)=^ + — 
I o 

V 



F2iv)-F2{-v)=jW + T^'^] 
F3(-q) + F3{-rj) =-^+2^1' 



9 1 2 



(A9) 



The number densities, energy densities, and fluxes for all other fermionic species (pairs and X-type neutrinos) can be quickly 
obtained repeating the procedure outlined above. 



APPENDIX b: conditions for silk damping, diffusive heating, and damping from lepton pressure loss 

In order to see if diffusive heating occurs at all for short wave-length compressible waves, the "Silk determinant" must be greater 
than zero i.e., 

fds\ fdU\ U + P fdYL\ 
«r — =7^ M.^'W^ >0. (Bl) 

In order to evaluate this expression, we need to specify an equation of state. Photons and leptons obey relativistic equations of 
state for bosons and fermions, respectively. However, the quantum mechanical behavior for the non-relativistic baryons varies with 
increasing depth. In the inner and outer core, the baryons are mildly degenerate, while near the neutrino sphere baryons behave as a 
classical ideal gas (P cx nT). 

It is convenient to define a fiducial Fermi energy for the non-relativistic baryons 
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such that the Fermi energy and and zero-point pressure for the neutrons and protons can be conveniently written as 

e„,, = Y^!^\F and Pfl = \ull^=lY^lp\eF (B3) 

2/3 

respectively. At early times in the mantle and near the neutrinosphere, ef ~ ^Pn MeV while ksT ~ 10-30MeV such that the 
contribution to the baryon pressure arising from degeneracy may be neglected at first approximation. Deep in the core below 



the mantle, the degeneracy pressure of the baryons is substantial such that e„,p 
15 - 20Mev (Pons et al. 1999). Their pressure is approximately given by 

2 /, 2- 



2/3 

F„,p 60 MeV is significantly larger than ksT 



3(0) 



1 + 



(B4) 



where we have kept the substantial finite temperature contribution for each baryonic species. The only species significantly con- 
tributing to the first term in eq. (Bl) are the neutrons and electrons where 
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Terms OiksT / pe) have been dropped for the sake of simplicity, which is why the proton contribution has been neglected. This is not 
a bad approximation during early times throughout the core where pg/kgT ~ 10-20 (Burrows & Lattimer 1986, Pons et al. 1999). 
The last term in eq. (Bl), which is largely responsible for diffusive heating of compressible perturbations is given by 
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Thus, to OiksT / pe) the Silk determinant is given by 

'ds\ 5YnUP^ r Stt^ 
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By putting in typical values for the core (F,, --0.7, F/, --0.34, e/r --60MeV, kgT ^ 20MeV, ^J-e-fJ-u, ^ 30MeV) the Silk determinant 
is approximately given by 



nT 



ds ' 
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(B9) 



which of course, is greater than zero. In order to realistically assess this result, we must note that the typical energy scales in this 
problem are given by the thermal energy per particle ksT ~ 20MeV and the Fermi energy Ef ^ 60MeV. Quantitatively, the value of 
the SiUc determinant given above is obtained by adding and subtracting characteristic energies of order ~ 20 -60 Me V. Therefore, 
the answer given by eq. (B9) should be taken with extreme caution and a more careful calculation is necessarily required. 

Even though the Silk determinant is positive, a given compressible wave is not necessarily driven. Unstable driving from the 
entropy perturbation must overcome the loss of degeneracy pressure due to lepton diffusion during compression. From the form of 
the pressure perturbation given by eq. (56) or equivalently from the asymptotic growth rate given eq. (83), SiUc driving or diffusive 
heating will occur if 
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Derivatives of the total pressure with respect to temperature and neutrino chemical potential are approximately given by 
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Driving from radiative diffusion therefore occurs if 
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which is marginally satisfied for ksT > 20 MeV and/or p < 10 gem 



APPENDIX C: CHEMICAL AND THERMAL EQUILIBRIUM 
CHEMICAL AND THERMAL SOURCE TERMS 

In order to calculate the conditions for chemical and thermal equilibrium, we must linearize the Boltzmann equation for all six 
(n,p,e~,e^,i^e,i'e) particle species. Electron (beta) capture is the inelastic radiative process that mediates chemical and thermal 
balance i.e., 

e~+p ^ n + Ve (Dl) 
e^ + n ^ p + i'e- (D2) 
In order to close our set of equations and for the sake of simplicity, we make five simplifying assumptions before doing so: 

• All species are represented, to lowest order, by a Fermi-Dirac distribution function characterized by an individual temperature 
and chemical potential. For example, in the case of protons 

fp= (D3) 
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• Thermal equilibrium for the background state, 

• Chemical equilibrium for the background state. 



Tp — Tyi — Tg — Tg+ — Ty^ — Ty^ — T. 
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Due to strong and electromagnetic interactions 



STe- = STe+ = STp = ST„ = STg 



and 



The last assumption is that the medium is sufficiently hot and optically thick such that Ee :^ Ei, >> {m„-mp)c^. 
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Now, consider the evolution of the partial fraction Y^- mediated by the reaction e + p n + Ve, governed by the Boltzmann 
equation 

(dY,-\ Gj {8l + 3gl) lire 2 f , 2 f , r r,2. ^2.,, ..w, ,n8^ 
V /capture {Ticf {2ttY {TicY J h^c^J Jo 



where £0 is the electron energy. At the same time, we must consider electron emission via neutrino capture via the reaction Ve + n- 
e~ + p. This gives us 
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After a quick algebraic manipulation, we may write the familiar expression 
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Note that the term in the exponential on the right hand side is equal to zero in equilibrium. 
Now we can consider perturbations. By applying energy conservation we have 
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If we neglect the final state blocking of the baryons, we are left with 
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This integral can easily be calculated by making use of the relation 
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The expression for electron evolution becomes 

cCTo f kT V fkT 



6Tg-ST^^ 



■y/kT 



[F5{Ve)-Fs{v.j] + 



^S^ip + S^e -S^„- 5^„^ _ fi^^, 5Tg- STjy^ 



kT 



kT 



[F4iVe)-F4{Vi^j] 



where 



^O^cro{gl + 3gl) 
2 TTCTO (fief 

- —rz — zzi 



F„(r])= / dx 



x" 



(D9) 



(DIO) 



(DID 



.(D12) 



(D13) 



(D14) 



(D15) 



e-^-'' + 1 

If we wish to follow the evolution of the total electron fraction Y^ = Yg- - Ye+ where the evolution of the e'^'s are mediated by the 
reactions e"*" + n p + z7(,, then the relevant expression becomes 
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and the expression for conservation of electron neutrino fraction, which is defined as = - i^,, reads 
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where we have defined 

D„ = F(ri)-F(-r,). (D18) 

The energy source term for the electrons is calculated in a very similar manner To be clear, the first law of thermodynamics for 
the electrons may be written as 

ut TT{nc) (2tt) {Tic) 



■ j d^Pp j d^Pn dEodnoE^pil-f,,)/, (l -/^J x ^1 - 



Upon perturbations, conservation of energy per unit volume for the the electrons is expressed by 



(D20) 



Note that the source term which drives electrons into thermal equilibrium is not exactly equal to the source term which drives the 
electrons into chemical equilibrium even though they are comparable in both appearance and magnitude. This is of importance 
when considering a spectrum of perturbations in the upper layer of a proto-neutron star. For a given dynamical timescale, the depth 
beneath the neutrinosphere at which the the fluctuations are in chemical equilibrium may differ, in principle, from the layer in which 
thermal equilibrium is achieved. In terms of neutrino bubble driving, whether or not the perturbations are either in chemical or 
thermal equilibrium or both, will determine whether or not driving due to changes in neutrino chemical potential (or the lepton flux) 
or driving due to changes in neutrino temperature (or the radiative flux) dominate the mechanics of the instability. 
The first law for the combined electron, positron, neutron, and proton gas is given by 
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where 
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Pg=P„ + Pp + P, and P, = P,-+Pe* 

C„{ri) = F„{v)+F„{-ij). (D22) 
For the electron type neutrino's, the first law of thermodynamics reads 
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where now Ui,^ = U^^, + U^^ and P^^ = P^^^ +Pv^,- Notice that we have not included the non-electronic type neutrinos in this discussion. 
The only mechanism which thermodynamic ally couples the /i- and r-type neutrinos to the other fluid and radiation components is 
inelastic scattering with electrons, in other words, the inverse Compton effect. This coupling is heavily blocked by the final state 
occupancy of the electrons which reduces the scattering cross-section by a factor of (1 - /,.). 

CHARACTERISTIC TIMESCALES FOR CHEMICAL AND THERMAL EQUILIBRIUM 

Timescales for both chemical and thermal equilibrium dictate whether or not the source terms represented in eqs. (D16) and (D21) 
vanish. That is, electron capture and emission dictate both chemical and thermal equilibrium simultaneously. 

In order to calculate the magnitude of the lepton source term in eq. (D16), it is convenient to write the left hand side as 



-iLo5Y, + 5\ ■VY, = -, 
Thus, the source of electrons will vanish if 
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and 



where we have defined 

_ ao / kT 
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he 

where Ca is the characteristic absorption opacity and «o is a characteristic relativistic particle number density. Note that the expression 
for th was obtained by rewriting the condition for chemical equilibrium as 

Sj^,±S^^^-S^^ (E5) 
ki kT T 

Likewise, in order to calculate the net effect of the energy source term in eq. (D21), it is useful to write the left hand side as 



-/(5C/, + (5vVf/, + /(C/,+P,)k.,5v = -/ ( ^r,+ (:^) 5t^,+ [^] 5p 

The heat source for the yipe^ gas will vanish if 
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and 
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where 

f/o = 4— (E9) 

serves as a characteristic energy density for relativistic particles. The equilibrition rate t^F.chem, given by eq. (E3), is identical to that 
derived by Burrows & Lattimer (1986, their eq. 16) except for the presence of a factor {dYe/ drie)T,,p in the denominator and the 
absence of a term that drives the system into statistical equilibrium in the numerator. Both discrepencies can be readily accounted for 
by realizing that our analysis is concerned with perturbations about an equilibrium already in statistical and thermal balance, while 
those authors calculate whether or not their equilibrium is in statistical and thermal balance with respect to the contraction timescale. 

In what follows, we approximate whether or not acoustic-type fluctuations near the neutrinosphere are in chemical and thermal 
equilibrium. We conclude that chemical and thermal balance between the radiation and gas is marginally satisfied at the neutri- 
nosphere. 

REQUIREMENTS FOR CHEMICAL AND THERMAL EQUILIBRIUM NEAR THE NEUTRINOSPHERE 

Interestingly, the chemical and energy source terms which lead to the characteristic frequencies tjy th and tJu.th vanish. This results 
from the mathematical fact that F„{ri) is an even function if n is odd and an odd function if n is even. From a more physical point of 
view, chemical and thermal balance simultaneously occur as long as eq. (E5) is equal to zero. 

Throughout the star, from the core to the neutrinosphere, the electrons can be considered fully degenerate such that r/^ >> 1. The 
thermal contribution to the electrons is more profound near the neutrinosphere, but there, rje is still > a few. Thus, the Fermi-Dirac 
integrals for the electrons Fnirje) can be approximated as (Pathria 1996) 

l(n + 2) n+l 

especially if n is large. Note that r(n) is the Euler Gamma function. Also, near the neutrinosphere, the chemical potential of the 
electron-type neutrinos fi^^ ^ a few times smaller than fig. As a result, we may write 
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Putting in numerical values typical for conditions at and around the neutrinosphere allows us to write the characteristic frequencies 
for chemical and thermal equilibrium in terms of Y^, T, and fi^, 
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where Yij-i, T4, and /^20 are the electron fraction in units of 1 /3, temperature in units 4 MeV, and electron chemical potential normal- 
ized to 20 MeV, respectively. In the above formulae, we made the approximation that the gas is primarily supported by the thermal 
pressure of the baryons. At the neutrinosphere, the sound crossing time ~ 0.1 ms, implying that chemical and thermal equilibrium is 
marginally satisfied for acoustic-type oscillations. 



If the atmosphere is dominated by magnetic pressure, the acoustic time for the fast magnetosonic wave is shorter than that of ordinary sound waves. Even if this 
were the case, the slow magnetosonic wave is the excitation species of interest since it always has larger growth rates than the fast wave (BS03). 



